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ABSTRACT 



A bi-quadratic isoparametric plate/shell bending finite element is developed 
to study the behavior of isotropic and laminated composite plates. The element 
is based on Mindlin-Reissner’s theory and the principle of virtual displacements. 
The element is implemented in a computer program. Results are presented and 
compared with analytical solutions to validate this element. Good agreement is 
observed for thin plates, while discrepancies are noted for thick plates. Effects of 
various integration schemes on the element performance are presented. Convergence 
studies for laminated composites for different fiber orientations are also discussed. 
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I. INTRODUCTION 



A. THE SCOPE OF THE THESIS 

Finite element analysis provides a general tool to solve problems in struc- 
tural mechanics. The methodology is applicable for static and dynamic response of 
structures and in predicting the elastic stability limits. 

The focus of the present study is to develop tools to analyze laminated com- 
posite plates and validate the model by comparing with known solutions. 

More specifically, the objectives of the present study are: 

a) to review some of the pertinent literature in the area of laminated 
composite plates. 

b) develop a finite element for the analysis of composite plates. 

c) develop consistent mass and load matrices. 

d) study the effect of thickness to characteristic length of the plate. 

e) study the effects of integration schemes. 

The outline of the remainder of this thesis is as follows: 

The basic formulation of the stiffness, mass and load matrix for the bending 
of flat plates using laminated composite materials are described in Chapter II. 

Chapter III addresses certain aspects of computational implementation. 

Chapter IV describes some test cases, example calculations and comparison 
with classical plate theory. 

Finally, Chapter V reflects experience gained and some suggestions for future 
research. 
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B. LITERATURE SURVEY 



The finite element method, [Ref. 1], may be described as a general discretiza- 
tion procedure of continuum problems, posed by mathematically defined statements 
with applications to several engineering analysis problems. 

A brief literature review pertaining to the analysis of plates/shells using finite 
element approximation is presented in the following paragraphs. A significant con- 
tribution to include shear is given by Mindlin [Ref. 2], while Hughes etal, [Ref. 3] 
adapt this theory to develop finite elements for the analysis of isotropic plates [Refs. 
3, 4]. Laminated plate theory based on the classical KirchofF hypothesis has been 
established by Reissner and Starsky [Ref. 5] and Whitney and Leissa [Ref. 6]. The 
effect of reduced integration in isoparametric elements was presented by Zienkiewicz 
etal [Ref. 7] and Hughes etal [Ref. 3]. 

The finite element method of analysis for the plate bending problem including 
shear deformation has been presented by Pryor and Barker [Ref. 8]. Mawenya 
describes formulations for multi-layer plates [Refs. 9, 10]. The higher order shear 
deformation theory of laminated composite plates was developed by Krishna Murty 
[Ref. 11], and Lo et.al. [Refs. 12, 13] present a higher order, three-dimensional 
theory. 

Burt [Ref. 14] presented a higher order theory and compared with Pagano’s 
elasticity-theory solution for the case of cylindrical bending and a symmetric cross- 
ply laminate consisting of three equal-thickness layers. Bending of simply supported 
thick rectangular plates was presented by Srinivas and Rao [Ref. 15]. Exact elas- 
ticity solutions for some particular plate bending problems have been obtained by 
Pagano [Refs. 16, 17, 18, 19] and Srinivas and Rao [Ref. 20]. 



2 



Application of classical shell theory, including transverse shear deformation is 
presented by Vinson and Chou [Ref. 21]. Naschie [Ref. 22] studied large deflec- 
tion behavior of orthotropic composite materials. Srinivas [Ref. 23] developed a 
refined approximate theory for the static and dynamic analysis of finite, laminated, 
composite, circular cylindrical shells with general boundary conditions. 

Plate theories, which include shear deformation has been developed by Whit- 
ney [Ref. 24] and Mau [Ref. 25]. The first such theory for laminated isotropic 
plates is due to Yang, Norris and Stravinsky [Ref. 26]. Reddy [Ref. 27] developed 
a higher order shear deformation theory of laminated composite plates. Reddy and 
Sandidge [Ref. 28] presented mixed finite element models of the classical and shear 
deformation theory. The effect of transverse shear deformation on bending of elastic 
symmetric laminated composite plate undergoing large deformation is presented by 
Gorji [Ref. 29]. 

Based on anisotropic elasticity, Hearmon [Ref. 30] and Lekhnitskii [Ref. 31] 
present general theories for laminates. The covariant form for the transformed lam- 
ina stiffnesses has been given by Tsai and Pagano [Ref. 33]. Gibert and Schneider 
[Ref. 34] directed their study in this direction. Noor and Mathers [Refs. 35, 36] pre- 
sented the effects of shear deformation and anisotropy on the response of laminated 
anisotropic plates. 

Nelson and Lorch [Ref. 40] compare the accuracy of various plate models to 
predict the behavior of laminated orthotropic plates. 
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II. THEORETICAL FORMULATION 



A. INTRODUCTION 

In this chapter, the derivation of plate finite elements based on Mindlin’s 
theory is described. The principle of virtual displacements is invoked to obtain 
equilibrium relations. 

B. THE PRINCIPLE OF VIRTUAL WORK 

In this section, we prove that total internal virtual work is equal to total ex- 
ternal virtual work and equivalence of this principle to the minimum total potential 
energy principle. 

In general, the total potential energy of a structural system is equal to the 
sum of strain energy and potential energy. 

n p = u + v (2.i) 

where, 

lip : Total Potential Energy of Structural System 

U : Total Strain Energy 

V : Total Potential of External Loads 

The total minimum potential energy requires that first variation of total po- 
tential energy be zero, or 

Slip = 0 (2.2) 

or, 

SU + SV = 0 (2.3) 
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in other words, 



6U = -6V (2.4) 

This may also be written in a different form, recognizing U as being the work 
done by internal forces and that work done by the external forces being equal to the 
negative of the total potential energy of the external loads. That is, 



sw int = 6W ext (2.5) 

which is a statement of the principle of virtual work, stating that, if the body is in 
equilibrium, the total virtual work done by the internal forces is equal to the total 
virtual work done by external forces for arbitrary, kinematically admissible virtual 
displacements. 

It may be noted that the form in Equation (2.4) is restricted to conservative 
loadings while the form in Equation (2.5) is applicable for any loading form. 

The total internal virtual work may be written as 

8W ini = / {a} T {6t}d(A) (2.6) 

Jvol 

where, 

{<t} : Vector of Stresses 

{£e} : Vector of Virtual Strains 

By using generalized Hooke’s law for material constitutive relations, stresses 
may be expressed as 



W = [< Q ){* } 

where the matrix [Q] contains the material stiffness coefficients. 



(2.7) 



If the thickness t is constant, then total internal virtual work takes the form, 
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&W iM = I !{ £ } t [c?j {6c}d(vol) (2.8) 

J A 

In order to derive the element matrices, the principle of virtual work, which is 
equally applicable to the element as well as to the total structure, is applied to the 
element. The virtual work is additive and results in the virtual work of the entire 
structure under consideration. 

The linear strain-displacement relations are given by 

{e} = [5]{u} (2.9) 

while the virtual strains are given by 

{6c} = [ B } {6u} (2.10) 

The operator matrix, [5], is dependent on the shape functions and their deriva- 
tives, and {u} and {6u} are vectors of displacements and virtual displacements, 
respectively of the element. 

On substituting Equation (2.9), (2.10) in Equation (2.8), we obtain, 

SW M = I t({Su} T [B) T )[Q)(\B]{Su})dA (2.11) 

J A 

or, rearranging 

8W int = {8u} T (t [ [B} T [Q][B]dA) (2.12) 

J A 

SW in , = {<u} T |A'l (2.13) 

where [K] = t j A [B T ][Q][B\dA is the element stiffness matrix. 
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In order to derive mass and load matrices, we start from the expression for 



external virtual work, 



SW ext = J T{6u}ds - J x{6u}d(vol) + {<5u} T {F} (2.14) 

where, T is the external force per unit length along the boundary of the element, 
x is the body force per unit volume (inertia force) 

{F} is the point loads applied at nodal points. 

On substituting for displacements in terms of nodal displacements using shape 

functions, we obtain, 



6W ext = { 6u} T (J[N]Tds ) - {^} T ( J tp[N}[N] T dA)[u) + {Su} t {F} (2.15) 

By invoking the principle of virtual displacements and noting that {<5u} is 
arbitrary, we get the equilibrium equations in the following form; 



[*-]{«} + [A/]{u} = [f] + [fr ( 2 . 16 ) 



where, 



[M]= f tp[N][N] T dA 

J A 



(2.17) 






[F ] = { \ 



F n 



(2.18) 



[F] d = f\N\Tis (2.19) 

It may be noted that the matrix [ M ] is the mass matrix, [F] is the vector of 
point loads and [F] d is the vector of consistent loads. 



7 



C. LAMINATE THEORY 



1. Introduction 

In the expression for [K\ matrix, there are three unknown matrices. In 
this section, we discuss about calculation of [Q] matrix. The matrix [Q] , relat- 
ing stresses and strains, consists of material stiffness coefficients. It reflects the 
properties of both fibers and matrix. 

The behavior of laminated plates is characterized by possible coupling 
between membrane action and bending action. 

Following discussions review some aspects of the laminate theory. 

2. Strain-displacement Relations 

As shown in Figure 2.1, a general strain field converts configuration 012 

to 0T2'. 

Linearized normal strains may be written as: 

du 

€x = fa 

and 

dv 

The shear strain is defined as the amount of change in a right angle. For 
small angles, 





Ixy — fil "h 02 

du dv 
dy ^ dx 



( 2 . 22 ) 

(2.23) 



Similar expressions may be derived for other strain components. 
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The foregoing strain-displacement relations can be summarized in a matrix- 



operator form. 



> 

£ X 








< e * 




7 xy 




lyz 




. 7 ZX j 





d_ 

dx 

0 



0 
a_ 

dy 

0 

d_ JL 
dy dx 



0 



o -2- 



d_ 

dz 

0 



1 
dx J 



3. Stress-Strain Relations 

The generalized Hooke’s Law takes the form, 



W = [<?]{«} 



where, the stress vector is given by 



(2.24) 



(2.25) 



{ G } — { G x Gy G 2 ~xy Gy Z G zx } 

and the strain vector is given by 



(2.26) 



{c} — {e* Cy c z 'Yxy ‘jyz Hzx} (2.27) 

The material stiffness matrix [Q] , contains nine independent coefficients. 

An orthotropic material has only five independent elastic coefficients —Ex,Ei, ^ 12 , ^ 21 , and G\i- 
More explicitly, Equation (2.25) may be written as 



> 














r ' 


°x 




Qll Q \2 Ql 3 


0 


0 


0 




€x 


°y 




Q2I Q22 Q23 


0 


0 


0 




e y 


G z 


► = 


Qsi Q32 Q33 


0 


0 


0 


< 


tz 


Try 




0 0 0 


Q 44 


0 


0 




7 xy 


T yz 




0 0 0 


0 


Q 55 


0 




lyz 


Tzx , 




1 

0 

0 

0 


0 


0 






. 7 ** > 



(2.28) 
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Figure 2.1: Displacement and Distortion of Differential Lengths dx and 
dy. 
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For planar orthrotropic material, the stress-strain relation reduces to: 



> 
















y 


< 7 x 




Q 11 


Q 12 


0 


0 


0 




c x 


a y 




Q 21 


Q 22 


0 


0 


0 


J 




T X y 




0 


0 


Q44 


0 


0 


\ 


7 * 1 / 


T yz 




0 


0 


0 


Q 55 


0 




7 yz 


L T ZX J 




. 0 


0 


0 


0 


Qes . 




, Izx J 



where, 



(2.29) 



Q 11 

Q 22 



E 1 

1 - i/ 12 V 2 \ 

E 2 

1 — ^12 ^21 



$44 



£*12 



Q55 — £*23 



<?66 — £*13 



(2.30) 

(2.31) 

(2.32) 

(2.33) 

(2.34) 



It may be noted that shears r y2 and t zx are obtained to account for 
transverse shears. 

4. Lamina of Composite Materials 

Composite structures are built of individual lamina, which are stacked 
into several number of layers to form a laminate. Each lamina consists of, typically, 
uniaxial fibers embedded in a matrix, such as a resin. In Figure 2.2, the principal 
material axes are labelled 1 and 2, that is, 1 -direction is parallel to the fibers direction 
and the 2-direction is normal to them. 

It may be noted that in each lamina there exists a state of plane stress. 
The state of stress is also shown in the Figure 2.2, in both 1-2 and x-y 
coordinate system. The computation stresses in different coordinate system follows 
the usual transformation rules [Ref. 21]. We follow the notation that <r 1 , rr 2 are 
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Figure 2.2: Lamina Coordinate System (2-Dimension) 
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normal stresses, <73,0-4 and <75 are the shear stress in 1-2 system. ei,e2,^3,£4 and e 5 
are the corresponding strains in 1-2 system. 



CS 

b b 

1 


= [T)cl 


H ;» 

b b 

1 


. . 




_ Tx v . 



where the transformation matrix is given by 



(2.35) 



m 2 


n 2 


2mn 


n 2 


m 2 


—2mn 


—mn 


mn 


(m 2 — n 2 ) 



[T\cl = 

The direction cosines of the unit normal are determined from 



(2.36) 



m = cos 9 



(2.37) 



n = sin 6 



(2.38) 



The subscript CL refers to two-dimensional case, that is the x-y plane 
only. Similarly, the strains are related in the two coordinate systems by 



lows: 



_ Cl 
e 2 


= r Ad 


€ x 

e y 


e 12 




- 7xy _ 



(2.39) 



By including transverse shears, we modify these transformations as fol- 



O' 1 




(7x 


(72 






(73 


= m 


7 X y 


O4 




T VZ 


- . 




. T zx . 



(2.40) 
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(2.41) 







^X 


£2 






^3 


= m 


7 xy 


e 4 




7 yz 






7 ZX 



where the transformation is given by, 



[T} = 



m 2 


n 2 


2mn 


0 


0 




n 2 


m 2 


— 2mn 


0 


0 




—mn 


mn 


m 2 — n 2 


0 


0 


(2. 


0 


0 


0 


m 


— n 




0 


0 


0 


n 


m 





The stresses and strains in x-y system may simply be obtained from 1-2 



system by inverting the matrix [T], 



and, 



with 



form, 



a x 




( \ 
o-i 


°y 


= PI" 1 < 


^2 


T X y 


03 


T V: 




C 74 


T Z x 




, . 



(2.43) 



€x 




ei 


ty 


= pt 1 


^2 


Ixy 


e 3 


lyz 




^4 


7 zx 




. ^5 



(2.44) 




m 2 


n 2 


—2mn 


0 


0 


n 2 


m 2 


2mn 


0 


0 


mn 


—mn 


m 2 — n 2 


0 


0 


0 


0 


0 


m 


n 


0 


0 


0 


— n 


m 



(2.45) 



The stress-strain relations, then, in x-y system assumes the following 
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(2.46) 
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cr x 




Qu 


Q 12 


Q 13 
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Q 21 


Q 22 


Q 23 
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* T xy 
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Q 31 


Q 32 


Q 33 
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Ix-y 


T V* 




0 
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Q 44 
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7yz 


. T 2X . 




0 


0 


0 


0 


Qss 




. 7z* . 



where, 



[<?] = (T)- 1 [0] pi 



(2.47) 



5. Laminate Analysis 

As discussed earlier, the laminae are stacked to obtain a laminate. In this 
section, the theory associated with the mechanics of laminates is described. 

Consider a laminate composed of N lamina. For the k th lamina of the 
laminate, Equation (2.46) can be written as follows: 



■V 

v x 




( 


°y 


• =I«1k' 


ty 


* T~x y 


Ifxy 


T yz 




lyz 


. Tzx . 


K 


. 7*i , 



(2.48) 



K 



where all matrices must have the subscript K due to orientation of the particular 
lamina with respect to the plate x-y coordinates. 

The functional form of the displacement for a plate are given by: 



u(x,y,z) = u 0 (x,y) + za(x,y) 


(2.49) 


v(x,y,z) = v a (x,y) + z0(x,y) 


(2.50) 


w(x,y) = w 0 (x, y) 


(2.51) 



where u 0 , v 0 and w Q are middle surface in plane displacements, a and /3 are related 
to the rotations. 
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Substituting Equation (2.49), (2.50) and (2.51) into the Equation (2.24) 



results in: 



— 



^ = 



t, - 



£xy — 



€yz — 






du 0 da 

dx +z dx 


(2.52) 


dv 0 d/3 

~d^ + Z d^ 


(2.53) 


0 


(2.54) 


1 dvo\ 

2 \ dy dx J 


(2.55) 


1 / 5 dw \ 

2 ( /J+ aw) 


(2.56) 


1 / dw\ 

2 ( a+ a7) 


(2.57) 



The mid surface strains are given by the following relations: 



dug 

dx 
dv 0 

dy 

1 ( du 0 dv, 



e yo ~ 



e *y° — 0 





(2.58) 




(2.59) 


5l>o\ 




-fa 


(2.60) 



The curvature terms, associated with transverse bending are written in 



the form, 



da 
dx 

d£ 
dy 

1 (da a/3 ^ 
xy 2 \dy dx , 



— 



Ky = 



(2.61) 

(2.62) 

(2.63) 



On substituting the strain-displacement relations into the stress-strain 
relations, we obtain stresses in terms of displacement components, 
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(2-64) 



<?x 




t Xo + ZK X 


°v 




ty a + ZKy 


< T xy 


II 

S' 


~ixy 0 H" ZK xy 


Tyx 




~tyz 




K 


. 7** 



For plate/shell type structures, we define stress resultants N x , N y , N xy , 
Q x , Q y , M x , M y , and M xy , as shown in Figure 2.3. For the k t h layer, we may write, 



r N x ' 




\ 

O x 


Ny 


ft/ 2 


Oy 


< N X y 


= / . 


0 X y 


Qx 


J-t/2 


Oyz 


Qy J 




. ° zx > 



dz 



(2.65) 



For a laminate composed of N-layers, the normal stress resultants are 



given by 




or, in terms of strains, we have 



( 2 . 66 ) 




The equation (2.67) may be written in the form 



(2.67) 



(JV) = [A] M + (B] M (2.68) 

with the elements of the matrices [>t] and [B] given by 

N 

Aij = E(QiMh-t k . i) (i,j = 1,2,3) (2.69) 

k=l 

= (hi = 1.2,3) (2.70) 

Z fc=l 
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The moment resultants are obtained for a K tft layer as follows: 



M x ) 


ft / 2 1 


'*• 1 




< M y : 


= A 


<7y > zdz 


(2.71) 


, M X y J 


J-t/2 


. a xy ) 



[M] = [B] [e„] + [£>] [a] (2.72) 

where D {j (O.j) (f K ~ 1 k- i) (2.73) 

6 jt=i 

D. [5], THE STRAIN-DISPLACEMENT MATRIX 

1. Introduction 

In this section, we describe how the matrix [R] may be calculated. 

The matrix [R], which relates the strains and displacements at nodal 
points of an element may be obtained in the form, 

{e} = [R]{u} (2.74) 

It may be observed that this matrix depends on the choice of the nodal 
degrees of freedom, the shape function and the form of strain-displacement relations. 

In the discussion that follows, we describe the nine noded bi-quadratic 
Lagrangian isoparametric element. There are five degrees of freedom associated with 
each node. 

2. Shape Functions and Their Derivatives 

A finite element is called an isoparametric element if the same interpola- 
tion functions define both the geometry and displacements. 

Figure 2.4 shows the element in the mapped space. For the in-plane 
deformation, the geometry may be interpolated 
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Figure 2.3: Positive Directions for Stress 
for a Lamina 



Resultants and 



Stress Couples 
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(2.75) 



{;)- 



[N]{x 1 y 1 x 2 y 2 ••• x 9 y 9 } 



and the displacements given by 



| = [N]{u lVl u 2 v 2 ■■■u 9 v 9 ) (2.76) 

where the matrix of shape functions is given by 

Nr 0 N 2 0 ••• N 9 0 

0 Nr 0 N 2 0 N 9 

We may also write the displacements in a summation notation as 




(2.77) 




9 



£ Nu x 


(2.78) 


9 


£ N,v, 

t = l 


(2.79) 



In the case of plate bending, three more degrees of freedom are added to 
the planar displacements. 

The transverse displacement w , 0 X and 9 y are the rotations of the normal 
to the undeformed middle surface in the x — z and y — z plane, respectively. 



9 

w = £./V,in,- (2.80) 

:= 1 

9 X = X>,-0 xt (2.81) 

1=1 

0y = X>i*vi (2-82) 

t'=i 

The shape functions and their derivatives, which are needed for the com- 
putations of strains, are given in Table 2.1. 
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Figure 2.4 : Lagrangian Plate Element 
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Figure 2.5: Degrees of Freedom at a Typical Node i. 

It may be noted that the rotational degrees of freedom are treated as 
independent quantities, following Mindlin theory [Refs. 37, 38, 39]. 

3. Jacobian Transformation Matrix 

The chain rules for differentiation of N(£,tj) with respect to £ and tj gives, 



ON, _ d^dx dNj dy 
d( dx di **" dy d£ 



(2. S3) 



dN i _0N i dx dNj dy 
dy dx dr] + dy dr] 

or, in matrix form, 



' dN, ' 




dx 


Oy 1 




0N { -j 


dt 




di 


di 




dx 


dN, 




dx 


dy 




dNi 


dr] J 




. Or] 


On . 




dy J 



The Cartesian derivatives are given by 



(2.84) 



(2.85) 
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Figure 2.6: Displacements in the x-z Plane 
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TABLE 2.1: Shape Function and Their Derivatives 
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( 2 . 86 ) 



dNi 1 




dx 


dy 


-1 


’ dNi ' 


dx 






dt 




ft 


dNi 




dx 


dy_ 




dNi 


- dy . 




dr] 


dy 




dr] 



The Jacobian matrix [J], which contains the derivates of the global coor- 
dinates with respect to the local coordinates is given by, 



[J] = 



dx 


dy 


w 


dt 


dx 


dy_ 


dr] 


dr] 



(2.87) 



By using isoparametric element concepts [Ref. 1], we may write the Ja- 
cobian elements as, 



dNi 

dt X ' 


(2.88) 


dNi 

dt Vi 


(2.89) 


dNi 
dr] X ' 


(2.90) 


dNi 

dv Si 


(2.91) 



The inverse of the Jacobian matrix can be expressed as 

J 22 



i-i 1 



M “ U I i 

with the determinant of the Jacobian given by 



—J 12 



(2.92) 



J I— J\\J22 ~ J2lJl2 



(2.93) 
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The Cartesian derivatives are, then, given by 



II 

H 

co ^ 


1 

" \J\ 


(*.?■ 




(2.94) 


dNi _ 


1 i 


( dNi 


dNA 


(2.95) 


dy 


\J |l 


[ 21 dt 


+ j "en) 



4. Formation of [ B } Matrix 

The equation (2.74) may be written as 



Ixy 



JL o 

dx 

0 f 

d_ % 
dy dx - 



U 



(2.96) 



Replacing the displacements in terms of nodal displacements, we get 



or, in matrix form 




(2.97) 



W = [Bi] {u} 



where [5,] is given by 



[B.] = 

For plate bending, we have, 



“ dNj 
dx 


0 


o 




dNi 




_ dy 


dx 



t x 




1 

o 

81® 

> O 

i 




w 


ty 


— 


0 i f 

L u dy dx j 


< 


u > 








V 

\ 



or, in terms of nodal displacements, 



(2.98) 



(2.99) 



( 2 . 100 ) 
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or, 



E. 



1 

KO 

1 


— 


_ Try _ 





0 

0 

0 



d£N. 

dx 

0 

dJ^N, 



0 

aft* 



dy dx 

The transverse shears are expressed by 





r \ 




Ui 1 


< 


Ori 




Oyi ) 



lyz 

Izx 



f 0 1 

l 1 o 



lyz 




*)zx 





dx 



ENiWi 

- E Ni$xi 

-ZNiOyi 



-ZNi 



[Bi] = 



GAUSSIAN QUADRATURE 
1. Introduction 
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ZN t 




0 


i.te elements 


as 




0 


-dty 

dx 


0 


-1 


0 


0 




dNj 


o 






A 


u 


dy 




dx 


dNj 

dy 


0 




■Ni 


dNj 

dx 


-Ni 


0 


j 




(2.101) 



( 2 . 102 ) 



(2.103) 



(2.104) 



In this section, we discuss aspects of numerical integration. Here we 
describe the Gaussian method to calculate the integrals, which is widely used in 
finite element work [Refs. 38, 39]. 

2. Summary of Gauss Quadrature 
To approximate an integral /, 



/ = (2.105) 

where is a function defined in that region. As shown in the Figure 2.7, we 
sample <f> at the midpoint of the interval and multiply by the length of the interval. 
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Generalization of Equation (2.105) leads to the formula 



i = £ 

i=l 



(2.106) 



or, 



I = + W 2 <f>(( 2 ) + ■■■ + W n <j>tf n ) (2.107) 

where £, is the location of the integration point i with respect to the origin, is 
the weighting factor for point i and n is the number of points at which <^(£) is to be 
evaluated. 

The sampling points and weights for Gaussian Quadrature, which is 
adapted in the present study, are listed in Table 2.2. 

In two dimensions, we may approximate I by 



I = '£'£ W ‘ W >'KU,Vi) (2.108) 

» i 

In Equation (2.6) and (2.7) each coefficient in the integrand matrix must 
be integrated as described above. 
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Figure 2.7: Gaussian Quadrature (£ Coordinate) using Three Sampling 
Points 
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TABLE 2. 



n 


•if. 


R, 


1 


0 0 


2 0 


*> 


0 577350:692 


1 0 


3 


0 774590- '692 
0 0 


0 5555555556 
0 8K88KKS889 


4 


(1 KM 1 16 
0 VWKUU36 


0 3478548451 
0 6521451549 


5 


0 Vi >61798459 
0 53M093|01 
0 0 


0 2.M.926K85I 
0 -)7h62Kf>7()5 
0.5(>KKKK88R9 


6 


0 9324095142 
0 6612**9^65 
0.2*8619 1 K61 


0.1713244924 
0 3607615730 
0 4679139346 


7 


0 949KP9123 
0 741531 1S56 
0 4058451^14 
0 0 


0 1294849662 
0 2797053915 
0 3818300505 
0 417959 IS37 


8 


0 V602KV8565 
0 79n>664774 
0 5255324099 
0 1K34UM25 


0 1012285363 
0 2223810345 
0 3 1 37066459 
0 3626837834 



2. Coefficients for Gaussian Quadrature 
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III. PROGRAM IMPLEMENTATION 



A. INTRODUCTION 

After having formed the element stiffness matrices, the global matrices are 
formed in a standard manner [Refs. 1, 38], resulting in, for static analysis, 

{F} = [K] M 

where {F} : External loads vector 

{K} : Stiffness matrix 

{u} : Nodal displacements vector 

If we know {F} and [/<'], {u} may be computed for a static problem. Typical 
steps for static analysis may be described as follows: 

Step 1: Input the material properties, plate coordinates, loads, the number of 
integration points, boundary conditions. 

Step 2: For the composite, input number of layers, elements and fiber orien- 
tations. 

Step 3: Determine element matrices in global coordinate system. 

Step 4: Assemble the element matrices. 

Step 5: Solve for displacements and stresses. 

The Figure 3.1 shows the flow chart of the program. 

B. SOLUTION PROCEDURE 

In order to obtain numerical solution, the element matrices were programmed 
and incorporated into an existing finite element program. The implementation in 
double precision was done on a 32-bit Apollo D-3000 series computer. 

Following steps describe the element subroutine. 
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BEGIN 





Figure 3.1: Flow Chart (Main Program) 
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Step 1: Material properties from main program are read. 

Step 2: Calculate [Q] according to number of layers. 

Step 3: Select integration point 2 x 2 or 3 x 2 or 3 x 3. 

Step 4: Establish shape function. 

Step 5: Determine the Jacobian matrix and [£]. 

Step 6: Formulation of [K\. 

Step 7: For each integration point, do steps 4, 5, and 6, and accumulate [K\. 
Step 8: Calculate [ K] for each element. 

Step 9: Return to main program. 

The flow chart for element level computation is given in Figure 3.2. 
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BEGIN 




Figure 3.2: Flow Chart (KQML9) 
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IV. NUMERICAL EXAMPLES 



A. INTRODUCTION 

This section discusses the validation of the element formulation and imple- 
mentation by means of selected numerical examples. The isotropic elastic plates are 
solved under different boundary conditions and loadings. 

This is followed by application to selected laminated composite plates to check 
the formulation for such applications. The effects of thickness and integration 
schemes also are investigated. 

B. SIMPLY SUPPORTED SQUARE ISOTROPIC PLATE 

1. Simply supported plate under concentrated load. 

Consider a rectangular isotropic plate under a central point load. The 
geometry and the boundary conditions are shown in the Figure 4.1. The structure 
is modeled using the double symmetry. 

The properties of material / are given by 
E = 10.92 xlO 6 (psi) 

v = 0.3 
L = 10 in 

P = 400 (lbs) 

The geometric boundary conditions of the simply supported plate are 
w = 0 on all edges. 

As shown in the Figure 4.2, results depicting maximum displacement 
obtained using 2x2 integration points, 3x3 integration points and ‘heterosis’ [Ref. 
4] elements are compared. As the thickness decreases, the element is more effective 
for integration schemes and for thick plates, the error is approximately 20%. It is 
possible that by increasing the number of elements, it may improve the solution. 
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Figure 4.1: A Rectangular Isotropic Plate under Concentrated Point 
Load 
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WlW ana i vs. No. of elements are shown in Figure 4.3. As the number of 
elements are increased, the solution converges to the analytically predicted values. 
It may be noted that reduced order of integration is more flexible, and approaches 
the analytical solution as an upper bound. 

2. Simply Supported Plate under Distributed Load 

Next we consider a rectangular plate under a distributed load. The ge- 
ometry and the boundary conditions are the same as in the previous example. The 
model of the structure, using the symmetry, is shown in Figure 4.4. 

The material properties for this example are given earlier as material I. 
The normalized maximum deflection is shown plotted against the number 
of elements in Figure 4.5. It may be observed that different integration schemes 
converge towards a lower bound. Figure 4.6 depicts the effectiveness of this element 
in predicting behavior of thick and thin plates. 

This element gives better predictions for thick plates than heterosis ele- 
ment, however, for thin plates, heterosis element appears to be better. 

C. CLAMPED-CLAMPED SQUARE ISOTROPIC PLATE 

Consider a rectangular plate clamped on all sides, under a distributed load. 
The geometry and boundary conditions are shown in Figure 4.7. The boundary 
conditions for this problem are implemented by prescribing all degrees of freedom 
to zero, on all four edges. The structure, again, is modeled using the symmetry. 

The material properties for this plate are same as for material I. 

The results showing the normalized maximum displacements are shown in 
Figure 4.8. The heterosis element results in about ten percent error while the 
Lagrangian element gives about twenty to twenty-three percent for thick plates. 
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Figure 4.2: IF/W ana j vs. Log (L/T) on Simply Supported Square Plate 
under Central Point Load. 
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Figure 4.3: IF/ W ana [ vs. No. of Elements on Simply Supported Square 
Plate under Central Point Load. 
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Figure 4.4 
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A Rectangular Isotropic Plate under Distributed Load 
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Figure 4.5: \VjW( ana i ) vs. No. of Elements on Simply Supported Plate 
Under Distributed Load 
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Figure 4.6: W/W {anal) vs. Log(L/T) on Simply Supported Plate Under 
Distributed Load 
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Figure 4.7: A Clamped-Clamped Rectangular Plate Under Distributed 
Load 
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However, for thin plates, the predictions are about the same using either elements. 
Figure 4.9 shows the convergence characteristics of the element. 

D. CANTILEVERED ISOTROPIC PLATE 

Next example considered is a cantilevered plate. The geometry and boundary 
conditions are shown in the Figure 4.10. The material properties, (material II), are 
given by 



E = lx 10 4 (k/in 2 ) (4-1) 

v = 0.3 (4.2) 

L = 10(m) (4.3) 

t = 0.2(in) (4-4) 



The result of this example is shown in the Figure 4.11. As the number of 
elements is increased, the results converge, for all integration schemes. 

E. SIMPLY SUPPORTED LAMINATED PLATE 
1. Graphite-epoxy 

Consider a rectangular composite plate under a sinusoidal load. Two 

types of construction are used. 

Case I : Graphite- Epoxy 0 o /90°/90 o /0° 

Case II : Graphite- Epoxy 0°/ — 60°/60°/0° 

Material properties are given by 



E 2 
@22 
E 2 



0.6 

0.5 



"12 
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Figure 4.8: W/W( ana i) vs. Log(L/T ) for Clamped-Clamped Plate under 
Distributed Load 
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Figure 4.9: Nondimensional Deflection vs. No. of Element for 

Clamped-Clamped Plate under Distributed Load 
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Figure 4.10: Cantilevered Plate 
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Figure 4.11: W/W [anal) vs. No. of Element for Cantilevered Plate 
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(4.8) 



El 

E 2 

The plate is subjected to a transverse sinusoidal load of intensity 

. 7TX . 

q = q 0 sin — sin 
Jj 

This problem is studied to compare finite element results with those of 
shear deformation theory [Ref. 24] and elasticity solution [Refs. 16, 17, 18]. The 
results are in Figure 4.12 for Case I. The deflection obtained from the finite element 
solutions agree very well with the exact elasticity solutions. The solutions agree 
well for thin plates, while a large discrepancy between the present predictions and 
Reference 24 is observed. 

The maximum deflection is plotted as the number of elements are in- 
creased in Figure 4.13. The convergence trend may be noted as the number of 
elements are increased. In Figure 4.14, the results for the second type of composite 
(0°/ _ 60°/60°/0°) is depicted. 

2. Glass-Epoxy 

This section considers rectangular composite plate under sinusoidal load. 
The two types of construction are Case III, 0°/90°/90°/0 and Case IV, 0° / — 
60°/60°/0°. The plates are square and simply supported. 

The material properties are given by 



E 12 
E 2 


= 0.6 


(4.10) 


G 22 

e 2 


= 0.5 


(4.11) 


V\2 


= 0.25 


(4.12) 


E i 

Eo 


= 3 


(4.13) 
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Figure 4.12: Central Deflection W max of 4-Layered Square Plate 

(Case I: Graphite-Epoxy) 
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Figure 4.13: W max vs. No. of Elements for 4-Layered Square Plate 

(Case I: Graphite-Epoxy) 
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ure 4.14: W max vs. No. of Elements for 4-Layered Square Plate 



(Case 11: Graphite-Epoxy) 
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The plate is subjected to a transverse sinusoidal load of intensity. 



q = q 0 sm — sin — (4.14) 

Figure 4.15 shows the maximum deflection of this element and of shear 
deformation theory [Ref. 24] and elasticity solution [Refs. 16, 17, 18] for analytical 
solution. Good agreement is obtained for thin plates while discrepancies exist for 
thick plates. Convergence characteristics are shown in Figures 4.16 and 4.17. 
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(Case III: Glass-Epoxy) 
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Figure 4.16: \V max vs. No. of Elements for 4-Layered Square Plate 

(Case III: Glass-Epoxy) 
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Figure 4.17: lF mox vs. No. of Elements for 4-Layered Square Plate 

(Case IV: Class-Epoxy) 
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V. CONCLUSIONS 



This study is directed towards understanding the linear static analysis of plates 
composed of both isotropic and laminated composites. 

The formulation is based on the principle of virtual work. A finite element 
formulation is presented as a model for the analysis of laminated anisotropic plate 
structures. Several numerical examples for the isotropic elastic plates are solved for 
different boundary conditions and loadings. 

The results show that bi-quadratic isoparametric lagrange plate bending ele- 
ment is effective for analysis of laminated anisotropic plates. Numerical solutions 
agree well with analytical solution. Further, it is observed that as the number of 
elements are increased, the convergence to analytical solution is assured. 

The element predicts good results for thin plates while large discrepancies exist 
for thick plates and calls for further investigation. In general, a 3 x 3 integration 
seems to predict the solutions well. 

More numerical experiments need to be done regarding selective integrations, 
and apply to a variety of layer configurations for composite plates. Another area is 
in including the nonlinear terms for buckling problems. The free-vibration analysis 
to predict the natural frequencies and mode shapes is an obvious extension. 
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